function FWD()
% Analyze a forward protocol
% This returns the modified control worksheet
% since some analysis results are placed in the sheet, and the data are all
% corrected in advance (this used to call prepare data...)

% n is flag for spike analysis - omit if n is 1, do if n is 2

% add return of ivresult structure 3/2/99. ivresult contains:
% ivresult.timestamp = time of latest analysis
% ivresult.iss is mean current
% ivresult.vss is mean voltage
% ivresult.iss_n is mean current excluding spikes
% ivresult.vss_n is mean voltage excluding traces with spikes
% ivresult.vpk_n is peak voltage excluding traces with spikes
%
% Modified 2/10/2000 : changed so that changing rates in the middle of the
% run does not cause subsequent traces to be incorrectly analyzed.
% i.e., all time pointers are recomputed for every trace, depending on the 
% current rate...
% Affected also: find_spikes, spike_sum, etc.
%

try
   sf = getmainselection;
   if(sf > 0) 
      pflag = getplotflag;
      QueMessage('Forward analysis', 1); % clear the que
      for i = 1:length(sf)
         FWD2(sf(i), pflag);
      end;
   end;
catch
   watchoff;
   QueMessage('Error in Forward Analysis routine', 1);
end;

function FWD2(sf, plot_flag)

global VOLTAGE CURRENT DFILE

global CONTROL 

QueMessage('FWD - Starting');
dat = [];
time = [];
ivresult = []; % initialize it

[DFILE, err] = analysis_setup(DFILE, sf);

if(err ~= 0)
   return;
end;

% abstract general information
spike_thresh=number_arg(CONTROL(sf).thresh);
QueMessage(sprintf('FWD - Spike Threshold = %7.2f mV', spike_thresh));
protocol=deblank(lower(CONTROL(sf).protocol));

[records,pts]=size(CURRENT);

% get windows 
ts1=number_arg(CONTROL(sf).durho);
ts2=number_arg(CONTROL(sf).durs1);
ts3=number_arg(CONTROL(sf).durs2);
ts4=number_arg(CONTROL(sf).durs3);

method_flag = 0; % set according the the protocol that was run
% designate times according to protocol
switch(protocol)
case {'forward'}
   tdel=ts1; % delay to first stimulus
   tsw1 = ts1; % Begin spike window 1 measurements
   tsw1e = tsw1+ts2; % End of spike window 1 measurements
   tsw2 = ts1+ts2+ts3; % Begin spike window 2 measurements
   tsl = tsw2; % latency zero point (spike window 2 start)
   tsw2e = tsw2+ts4; % end of spike window 2
case {'forward2'}
   tdel=ts1; % delay to first stimulus
   tsw1 = ts1; % Begin spike window 1 measurements
   tsw1e = tsw1+ts2; % End of spike window 1 measurements
   tsw2 = ts1+ts2+ts3; % Begin spike window 2 measurements
   tsl = tsw2; % latency zero point (spike window 2 start)
   tsw2e = tsw2+ts4; % end of spike window 2
   method_flag = 1;
   mseq = seq_parse(CONTROL(sf).sequence);
otherwise
   QueMessage(sprintf('FWD - Error: protocol %s not recoginized - Exiting', protocol));
   return;
end

% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
RATES = DFILE.rate .* DFILE.nr_channel / 1000; % array of sampling rates, convert to msec

% spike window times
sw1b = floor(tsw1./RATES);
sw1e = floor(tsw1e./RATES);
sw2b = floor(tsw2./RATES);
sw2e = floor(tsw2e./RATES);

% first smooth the voltage out a bit
for i = 1:records
	fsamp = 1000/RATES(i); % get sampling frequency
	fco = 5000;		% cutoff frequency in Hz
	wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
	if(wco < 1) % if wco is > 1 then this is not a filter!
	   [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
	   vsmo(i,:) = DigitalFilt(b, a, VOLTAGE(i,:)')'; % filter all the traces...
	else
	   vsmo(i,:) = VOLTAGE(i,:);
	end
end


% get spike times
QueMessage('FWD - Finding spikes Window 1');
[fsl_w1, fisi_w1, nr_w1, st_w1]=find_spikes(DFILE, tsw1, tsw1e, spike_thresh);
% disp(sprintf('Rec: found %4d spikes', sum(nr_spikes)))

QueMessage('FWD - Finding spikes Window 2');
[fsl_w2, fisi_w2, nr_w2, st_w2]=find_spikes(DFILE, tsw2, tsw2e, spike_thresh);
if(isempty(st_w2))
   QueMessage('FWD - Failed to find spikes in test window', 1);
   tsw2
   tsw2e
   return;
end;
% disp(sprintf('Rec: found %4d spikes', sum(nr_spikes)))
fsl_w2x = fsl_w2(~isnan(fsl_w2));
fisi_w2x = fisi_w2(~isnan(fsl_w2)); % strip out the NaN's: but only those where first isi is not present either
												% to match st_w2 length
                                    
QueMessage('FWD - Finding ALL spikes');
[fsl_wa, fisi_wa, nr_wa, st_wa]=find_spikes(DFILE, 0, 1500, spike_thresh);
                                    
% Compute mean isi for traces in which both first and second traces gave more than one spike
isi_m1 = [];
isi_m2 = [];
isi_nc = [];
for i = 1:records
   isi_m1(i) = NaN;
   isi_m2(i) = NaN;
   k1 = find([st_w1.source] == i);
   k2 = find([st_w2.source] == i);
   if(~isempty(k1) & ~isempty(k2)) % in record i, both trials had spikes...
      if(nr_w1(i) >= 2 & nr_w2(i) >= 2)
         isi_m1(i) = mean(diff(st_w1(k1).latency)); % mean isi first trial
         isi_m2(i) = mean(diff(st_w2(k2).latency)); % mean isi second trial
      end;
      if(nr_w1(i) == 0 & nr_w2(i) >= 2)
         isi_nc(i) = mean(diff(st_w2(k2).latency)); % mean isi for trials with no first response
      end;
   end;
end;

QueMessage('FWD - Analysis complete');

% plot if figure is set
if plot_flag>=0
   h = findobj('Tag', 'FWD'); % check for pre-existing window
   if(isempty(h)) % if none, make one
      h = figure('Tag', 'FWD', 'Name', 'Forward Protocol Analysis', 'NumberTitle', 'off');
   end
   figure(h); % otherwise, select it
   clf; % always clear the window...
   fsize = 0.2; % in cm...
   msize = 3;
   
   set(gca, 'FontUnits', 'centimeters', 'FontSize', fsize);
   
	subplot('Position', [0.1 0.75 0.35 0.20]);   % scatter plot of spikes vs trials
   for i = 1:length(st_wa)
      u = length(st_wa(i).latency); % how many spikes?
      trial(i).rep = [zeros(u,1)+i]';
   end;
   plot([st_wa.latency], [trial.rep], 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 1);
   label_ax ('Time (ms)', 'Trial #', fsize);


% scatter plot, fsl vs trials
	subplot('Position', [0.1, 0.50, 0.35, 0.20]);   
   plot([st_w2.source], fsl_w2x, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', msize/2);
   grid;
   label_ax('Trial', 'FSL (ms)', fsize);
   xm1 = xlim;
   
   % scatter plot, fisi vs trials
   subplot('Position', [0.1, 0.25, 0.35, 0.20]);
   plot([st_w2.source], fisi_w2x, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', msize/2);
	label_ax('Trial','FISI (ms)', fsize);
  xlim(xm1); % match axes limits...
   
   subplot('Position', [0.1, 0.03, 0.35, 0.20]);
   bar([1:length(nr_w1)], nr_w1, 1);
   hfp = findobj(gca,'Type','patch');
   set(hfp,'FaceColor',[0.5 0.5 0.5],'EdgeColor','w', 'LineWidth', 1.5)
   
   xlim(xm1); % match axis limits
   label_ax('Trial', '# of Conditioning Spikes', fsize);
   
   subplot('Position', [0.55, 0.75, 0.35, 0.20]);
   fsl_h = [];
   if(length(fsl_w2) > 10) % determine whether there are single or mutiple trials...
      skipf = 10;
   else
      skipf = 1;
   end;
   tn = 1:skipf:length(fsl_w2);
   for i = 1:length(fsl_w2)/skipf
      fsl_h(i) = mean(fsl_w2(1+(i-1)*skipf:i*skipf));
   end;
   plot(tn(1:length(fsl_h)), fsl_h, 'k^', 'MarkerFaceColor', 'k', 'MarkerSize', msize);
   label_ax('Trial Group', 'Mean FSL (ms)', fsize);
   
   subplot('Position', [0.55, 0.50, 0.35, 0.20]);
   fisi_h = [];
   for i = 1:length(fisi_w2)/skipf
      fisi_h(i) = mean(fisi_w2(1+(i-1)*skipf:i*skipf));
   end;
   plot(tn(1:length(fisi_h)), fisi_h, 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', msize);
   label_ax('Trial Group', 'Mean FISI (ms)', fsize);
   
   
   subplot('Position', [0.55, 0.25, 0.35, 0.2]);
   plot(isi_m1, isi_m2, 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', msize);
   hold on;
   xp = zeros(length(isi_nc), 1)+0.5;
   plot(xp, isi_nc, 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', msize);   
   hold off;
%   xlim([0 max(isi_m1)]);
%  ylim([0 max([isi_m2 isi_nc])]);
   label_ax('Mean Conditioning ISI (ms)', 'Mean Test ISI (ms)', fsize);
   
   % plot textual information abstracted from analysis...
   subplot('Position',[0.55,0.03,0.35,0.18])
   r_uncomp = CONTROL(sf).access;
   axis([0,1,0,1])
   axis('off')
   text(0,.9,sprintf('%-12s R[%d:%d] %-8s (analyzed on: %s)',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol, date));
   text(0,.8,sprintf('%s', DFILE.comment));
   text(0,.7,sprintf('%s', CONTROL(sf).Commentary));
   orient landscape
   
   % control printing and closing of window for automatic runs
   % plot_flag = 1 creates plot and leaves it up
   % plot_flag = 2 creates plot but closes it when done
   % plot_flag = 2 creates plot and prints it and then closes it
   if (plot_flag == 2)
      print -dljet3;
   end
   if plot_flag >= 2  
      close
   end
end

function label_ax(xl, yl, size)
% label axes and automatic small size
xlabel(xl, 'FontUnits', 'centimeters', 'FontSize', size);
ylabel(yl, 'FontUnits', 'centimeters', 'FontSize', size);
set(gca, 'FontUnits', 'centimeters', 'FontSize', size);
return;